library(tdaTS) #our package
library(plotly) #3D plotting
library(fields) #image.plot functionThe simple 3D plot function using function plot_ly
function from the plotly package.
myplot_3D <- function(data)
{
minD <- min(data$x, data$y)
maxD <- max(data$x, data$y)
M <- max(minD, maxD)
resScale <- rbind(data,c(0,minD,minD),c(0,minD,maxD),c(0, maxD,minD),c(0,maxD,maxD)) ### add 4 points in the corners for scaling x and y axes
plot_ly(resScale, x = resScale$x, y = resScale$y, z = resScale$t,
type="scatter3d", mode = "markers", color = resScale$t)
}We propose, as a starting point, the following four data generators. You give an example for each of the generators.
We close a segment into a circle, maintaining the length of the curve (2 pi)
data <- segment_to_circle(n = 1000,
change = 0.5,
time_sampling = "unif")
# an example with additional noise
data$x <- data$x + rnorm(n= 1000, mean = 0, sd = 0.2)
data$y <- data$y + rnorm(n = 1000, mean = 0, sd = 0.2)
######## plot the data
myplot_3D(data)The circle is smashed into a segment equal to the diameter of the initial circle.
data <- circle_extinction(n = 1000,
change = 0.5,
time_sampling = "unif")
myplot_3D(data)we split a segment at the changepoint in two at the center point with
regular speed with final gap = 0.1 here.
data <- segment_to_segments(n = 1000,
change = 0.5,
gap = 0.1)
myplot_3D(data)An example with no change. We deform the shape in a ellipse, both rotationg and shifting it.
data <- circle_move_distortion(n = 10000,
rotation = 2,
X_rate = 0.5,
Y_rate = 3)
myplot_3D(data)In most example, there is the possibility to use parameters
X_rate and/or Y_rate to extend or contract the
shape and then to change the density of points by level (not their
number)
data <- segment_to_circle(n = 1000,
change = 0.5,
X_rate = 0.5,
Y_rate = 3)
myplot_3D(data)In following sections, we will filter the Diagrams using 3 possible thresholds:
birth <- 0.15
death <- 0.17
diag <- 0.1The birth threshold is a vertical threshold
The death threshold is a horizontal threshold
The diag threshold is a diagonal threshold
The objective would be to remove the noise part of each of the persistence diagrams.
We still need to find the right approach to filter the noise (considering filtering by ratio death/birth maybe?)
In 3 steps
result with no noise
result with additional noise
segment comparison with matrix of cross Wasserstein distance using dimension 1 elements
set.seed(7)
data <- segment_to_circle(n = 1000,
change = 0.5,
time_sampling = "unif")
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)data$x <- data$x + rnorm(n= 1000, mean = 0, sd = 0.1)
data$y <- data$y + rnorm(n = 1000, mean = 0, sd = 0.1)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)nb_levels <- 20
all_PD <- Generate_All_Persistence_Diagrams(data,
nb_levels = nb_levels)
################
#### thresholding
################
all_filter <- remove_noisy_points(all_PD,
birth = birth,
death = death,
diagonal = diag,
infinity = FALSE)
################
#### distances: Wasserstein
################
mydim <- 1 #choice for dimension in PD
t(all_filter[,1:2]) %>% as.matrix()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## time 0 4 6 6 9 10 11 13 14 15 16 17 19
## dimension 0 0 0 0 1 1 1 1 1 1 1 1 1
D_wasserstein_2by2 <- distances_Persistence_Diagrams(all_filter, type = "2by2", dimension = mydim)
par(mfrow = c(1,1))
D_wasserstein_2by2 <- t(D_wasserstein_2by2)
image.plot(D_wasserstein_2by2[,nrow(D_wasserstein_2by2):1], axes = F, col = grey(seq(0,1, length = 256)))
axis(1, at=seq(0,1,length.out = nb_levels), labels= 1:nb_levels)
axis(2, at=seq(0,1,length.out = nb_levels), labels= nb_levels:1)
segments(0,0.5,1,0.5, col = 2, lwd = 4)
segments(0.5,0,0.5,1, col = 2, lwd = 4)set.seed(7)
data <- circle_extinction(n = 2000,
change = 0.5,
time_sampling = "unif")
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)data <- segment_to_segments(n = 1000,
change = 0.5,
gap = 0.5)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)data <- circle_move_distortion(n = 10000,
rotation = 2,
X_rate = 5,
Y_rate = 1, time_sampling = "discrete",
nb_levels = 20)
data$x <- data$x + rnorm(n= 10000, mean = 0, sd = 0.1)
data$y <- data$y + rnorm(n = 10000, mean = 0, sd = 0.1)
Plot_All_Persistence_Diagrams(data,
birth = birth,
death = death,
diagonal = diag,
nb_levels = 20)nb_levels <- 20
all_PD <- Generate_All_Persistence_Diagrams(data,
nb_levels = nb_levels)
################
#### thresholding
################
all_filter <- remove_noisy_points(all_PD,
birth = birth,
death = death,
diagonal = diag,
infinity = FALSE)
################
#### distances: Wasserstein
################
mydim <- 1 #choice for dimension in PD
t(all_filter[,1:2]) %>% as.matrix()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## time 0 1 2 3 4 5 6 7 8 9 10 11 12
## dimension 1 1 1 1 1 1 1 1 1 1 1 1 1
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## time 13 14 15 16 17 18 19
## dimension 1 1 1 1 1 1 1
D_wasserstein_2by2 <- distances_Persistence_Diagrams(all_filter, type = "2by2", dimension = mydim)
par(mfrow = c(1,1))
D_wasserstein_2by2 <- t(D_wasserstein_2by2)
image.plot(D_wasserstein_2by2[,nrow(D_wasserstein_2by2):1], axes = F, col = grey(seq(0,1, length = 256)))
axis(1, at=seq(0,1,length.out = nb_levels), labels= 1:nb_levels)
axis(2, at=seq(0,1,length.out = nb_levels), labels= nb_levels:1)
segments(0,0.5,1,0.5, col = 2, lwd = 4)
segments(0.5,0,0.5,1, col = 2, lwd = 4)Need to think…
Need to look at the litterature
Threshold diagonal ? death over birth ratio ?
library(TDA)
circle <- function(n = 1000, center = c(0,0), radius = 1, noise = 0.5)
{
pos <- runif(n, min = 0, max = 2*pi)
x <- center[1] + radius*cos(pos) + rnorm(n, mean = 0, sd = noise)
y <- center[2] + radius*sin(pos) + rnorm(n, mean = 0, sd = noise)
res <- matrix(c(x,y), nrow = n, byrow = FALSE)
return(res)
}
data <- circle()
DiagAlphaCmplx <- alphaComplexDiag(X = data ,
library = c("GUDHI", "Dionysus"),
location = TRUE,
printProgress = FALSE)
plot(DiagAlphaCmplx$diagram, col = 1 + cumsum(DiagAlphaCmplx$diagram[, 1]))plot(data, col = 1,xaxt="n", yaxt="n",xlab="", ylab="", asp = 1)
one <- which(DiagAlphaCmplx[["diagram"]][, 1] == 1)
for (i in seq(along = one))
{
for (j in seq_len(dim(DiagAlphaCmplx[["cycleLocation"]][[one[i]]])[1]))
{
lines(DiagAlphaCmplx[["cycleLocation"]][[one[i]]][j, , ] + rnorm(n = 4,sd = 0.0), pch = 19, cex = 1, col = i + 1)
}
}